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I. Introduction 

The linear bicharacteristic scheme (LBS) was originally developed to improve unsteady solutions 
in computational acoustics and aeroacoustics [l]-[7]. The LBS has previously been extended to 
treat lossy materials for one-dimensional problems [8], It is a classical leapfrog algorithm, but is 
combined with upwind bias in the spatial derivatives. This approach preserves the time-reversibility 
of the leapfrog algorithm, which results in no dissipation, and it permits more flexibility by the ability 
to adopt a characteristic based method. The use of characteristic variables allows the LBS to include 
the Perfectly Matched Layer boundary condition with no added storage or complexity. The LBS 
offers a central storage approach with lower dispersion than the Yee algorithm, plus it generalizes 
much easier to nonuniform grids. It has previously been applied to two and three-dimensional free- 
space electromagnetic propagation and scattering problems [3], [6], [7]. This paper extends the LBS 
to the two-dimensional case. Results are presented for point source radiation problems, and the 
LDTD algorithm is chosen as a convenient reference for comparison. 


II. Two-dimensional Implementation 

Maxwell’s equations in bicharacteristic form for linear, homogeneous and lossy media in the two- 
dimensional TM case (taking d/dz = 0) are 
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where c = 1 /^/ye, and the electric and magnetic conductivities are given by a and a*, respectively. 
The characteristic variables are defined as P = D z — / H y , Q = l) z ■ ^ H y , R = l) z ■ H x , 
S = D~ — i H r to represent the +:/• and ±y propagating solutions, respectively. Equations (l)-(2) 
can be rewritten more concisely as 
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where a = a/e + a* / y, and b = a fe — a* / y. To develop the discretized algorithm for a two- 
dimensional system, the stencils of Ligure 1 are proposed for the 2D LBS. The upwind bias nature 
of these stencils is clearly evident. We discretize time and space as t = nAt, x = iAx, y = j Ay, 
and we also know that H x = c{R — S)/2 and H y = c{Q — P)/ 2. Using the stencils shown in 
Ligure 1, the resulting finite difference equations for (3)-(6) are 
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Fig. 1. Two-dimensional upwind leapfrog computational stencils for (a) -\-x; (b) —x; (c) -\-y and (d) —y propagating 
characteristics. 
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where PJd denotes the value for P at grid point (i , j ) and time level n and i/ r = cAt/Ax, v y = 
(Al/ Ay are the x and y Courant numbers. Note that the differences are taken with respect to 
the cell center, i.e. the coordinate (i,j) is located at the center of the cell. The procedure for the 
TE polarization is similar, except the characteristic variables are defined by P = D y • / II z , 
Q = D y — \H Z , R = D r — \H Z , and S = D x + \H Z . The update equations are almost 
identical, differing only in signs of the cross derivative terms. The 2D LBS includes the PML 
boundary condition in its formulation without any additional storage or complexity. Starting with 
the bicharacteristic Maxwell’s equations in free space [(l)-(2) with a = a* = 0], it can be shown 
using a complex coordinate transformation PML approach that the resulting equations are identical 
to (3)-(6). A full derivation of these equations is outlined in [9], Thus, the LBS automatically 
implements the PML when combined with the standard PML conductivity profiles: linear, quadratic 
or geometric [10], 


III. Lourier Analysis 


Various Lourier analyses of the two-dimensional LBS 
have already been completed [5], [6], [7]; therefore, 
only the important results and conclusions from these 
previous analyses will be reviewed in this paper. Most 
of the information presented is summarized from [5], 
The stability condition for the 2D LBS is v :r . v y < 
1/2. Although this stability limit is more restrictive 
than the standard 2D LDTD method, it is not partic- 
ularly troublesome because many LDTD simulations 
use a Courant number of 1/2 for improved accuracy. 
Lor a Courant number v = 0.4 and propagation angle 
of 45°, the numerical dispersion decreases smoothly 
with increasing grid resolution as shown in Ligure 2. 
Lrom this figure, we see that the LBS has approxi- 
mately 1/2 the phase error as LDTD. Generally, the 
dispersion error for the LBS grows as v — \ 0. When v = 
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Fig. 2. Phase speed error versus grid resolution N 
for FDTD method and LBS with v = 0.4 and 
a = 45°. 

1 /2, numerical dispersion is zero along the 




coordinate axes and is maximum at 45° as shown in Figure 3 for a grid resolution N = 10 cells/A. 
The optimal Courant number for the LBS is 1/2. This Courant number offers much lower dispersion 
for most all propagation angles except those near a 45° vector. For v < 1/2, numerical dispersion 
decreases as both grid resolution and Courant number are increased. Typically, LBS dispersion is at 
least 1/2 that of FDTD, and can be much lower in many instances. 


IV. Numerical Results 



Fig. 


3. Phase speed error versus grid propagation 
angle a for FDTD method and LBS with = 
1/2 and N = 10. 


To demonstrate the 2D LBS, we consider various 
canonical problems. First, we inject an incoming plane 
wave on the outer boundaries using the LBS, and let the 
algorithm propagate the signal through the grid using a 
total field formulation. This is done by specifying the 
incoming characteristic variable (P,Q, R or S ) on the 
appropriate outer boundary. For example, on the left x 
boundary, P is specified for all j coordinates at i = 1. 

We use a 7 1 x 7 1 free space grid, with a Ax = Ay = 1 
cm, which has a time step of At = 16.67 ps and 
the incident wave is a Gaussian pulse with FWHM 
of 35 time steps (or 0.58 ns). We specify the inci- 
dence angle as 180°, and the electric field after 160 
time steps is shown in Figure 4. Similar results can 
be obtained with other incidence angles. It is clear 
that the LBS easily allows specification of incoming plane waves in its fundamental algorithm. 

Next we move on to radia- 
tion from a point source in free 
space. This problem demon- 
strates that the algorithm can 
easily treat spherical waves and 
it also tests the PML bound- 
ary condition. Two concurrent 
grids are used in this problem, 
each having a cell size of 1 mm. 
The first is a small test grid of 
size 101 x 101 cells with an ad- 
ditional 10 cell PML boundary 
condition. This grid is centered 
within a large 501 x 501 grid, 
and the point source is located at 
the center of both computational grids. The time step is 3.3 ps, and an electric field point source 
is located at the center of both grids and the total number of time steps is truncated at 512, to al- 
low no reflection from the large grid outer boundaries to reach the field sampling points. The inner 
grid is terminated with PML for both FDTD and the LBS, and the large grid is terminated with a 
second-order Liao boundary condition for FDTD and a characteristic based boundary condition for 
the LBS. The electric field is sampled at the same two locations in both grids, which are located 30 
cells in the +x direction from the point source and then ±30 cells in the y direction in the smaller 
grid. The point source is located in the smaller grid at grid point (61, 61) and the two sample points 
are (61,91) and (61,31). Figure 5 shows the electric field at the upper sample point (61,91) in 
the smaller test grid (with PML) for point source radiation in free space. Note the agreement is 
excellent, and there are no reflections from the outer boundary due to the PML boundary condition. 
Similar results were observed at the lower sample point (61, 31). Next, we simulate the same point 
source radiation problem, but turn off the PML boundary condition and instead use the character- 
istic boundary condition. Figure 6 shows the results using for the LBS with and without the PML 
boundary condition. Note the reflections from the outer boundary are clearly visible for the no PML 
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Fig. 4. Propagating plane wave injected on outer grid boundaries at 180° inci- 
dence. 


case. 




V. Conclusions 


This paper has extended the Linear Bicharacteristic Scheme for computational electromagnetics 
to the two-dimensional case. It was demonstrated that the LBS has several distinct advantages over 
conventional FDTD algorithms. First, the LBS typically has much lower dispersion than the FDTD 
method. The LBS can also be made to have zero dispersion error in certain instances. Second, 
the LBS provides the PML boundary condition as part of the fundamental scheme, and does not 
require additional storage or computational complexity for the PML. Finally, the upwind biasing 
provides a more flexible generalization to unstructured grids, and for subgrid implementations. The 
two-dimensional LBS can easily treat dielectric and magnetic materials, and a dielectric surface 
boundary condition was implemented in a separate report [9], However, validation of the 2D LBS 
for dielectric and magnetic materials was not explored in the present work and should be the subject 
of future articles. The LBS appears to be a very promising alternative to a conventional FDTD 
algorithm for many applications. Higher order schemes are available for the 2D case [5]. Additional 
research and development is required to explore the full potential of this new technique. 
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Figure 5: Sampled electric field due to radiation from 

a point source in free space. 
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Figure 6: Sampled electric field due to radiation from 

a point source in free space for LBS with and without 
PML boundary condition. 





